{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import string\n",
    "#https://stackoverflow.com/questions/19726663/how-to-save-the-pandas-dataframe-series-data-as-a-figure\n",
    "import six\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbpath = \"../../fits/\"\n",
    "productpath = \"../../postfit_derivatives/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              roi quantile         R0  car (week 0)  car (week 1)  \\\n",
      "0       AA_Global     mean   3.494587      0.100231      0.100732   \n",
      "1       AA_Global      std   1.541273      0.034553      0.052667   \n",
      "2     Afghanistan    0.025   1.468418      0.016030      0.012981   \n",
      "3     Afghanistan     0.25   2.446412      0.038256      0.037566   \n",
      "4     Afghanistan      0.5   3.436230      0.064816      0.068497   \n",
      "...           ...      ...        ...           ...           ...   \n",
      "1188     Zimbabwe      0.5   3.479923      0.162644      0.220131   \n",
      "1189     Zimbabwe     0.75   5.524083      0.469393      0.579358   \n",
      "1190     Zimbabwe    0.975  12.009383      1.292761      1.383765   \n",
      "1191     Zimbabwe     mean   4.307799      0.293625      0.370387   \n",
      "1192     Zimbabwe      std   2.969391      3.959783      2.956896   \n",
      "\n",
      "      car (week 2)  car (week 3)  car (week 4)  extra_std  ifr (week 0)  ...  \\\n",
      "0         0.103493      0.105365      0.096048   1.173057      0.002379  ...   \n",
      "1         0.065951      0.076988      0.079714   0.420742      0.002232  ...   \n",
      "2         0.012331      0.012318      0.012676   1.231066      0.001199  ...   \n",
      "3         0.038226      0.039458      0.041217   1.431145      0.003156  ...   \n",
      "4         0.071853      0.074996      0.078477   1.545750      0.005573  ...   \n",
      "...            ...           ...           ...        ...           ...  ...   \n",
      "1188      0.248206           NaN           NaN   2.646293      0.034050  ...   \n",
      "1189      0.613056           NaN           NaN   2.996698      0.088881  ...   \n",
      "1190      1.390096           NaN           NaN   3.795204      0.389523  ...   \n",
      "1191      0.360313           NaN           NaN   2.688589      0.078405  ...   \n",
      "1192      6.643321           NaN           NaN   0.505418      0.470311  ...   \n",
      "\n",
      "      car (week 10)  car (week 11)  car (week 7)  car (week 8)  car (week 9)  \\\n",
      "0          0.104464       0.084199      0.057180      0.065890      0.070085   \n",
      "1          0.099795       0.083278      0.048009      0.060984      0.078586   \n",
      "2               NaN            NaN           NaN           NaN           NaN   \n",
      "3               NaN            NaN           NaN           NaN           NaN   \n",
      "4               NaN            NaN           NaN           NaN           NaN   \n",
      "...             ...            ...           ...           ...           ...   \n",
      "1188            NaN            NaN           NaN           NaN           NaN   \n",
      "1189            NaN            NaN           NaN           NaN           NaN   \n",
      "1190            NaN            NaN           NaN           NaN           NaN   \n",
      "1191            NaN            NaN           NaN           NaN           NaN   \n",
      "1192            NaN            NaN           NaN           NaN           NaN   \n",
      "\n",
      "      ifr (week 10)  ifr (week 11)  ifr (week 7)  ifr (week 8)  ifr (week 9)  \n",
      "0          0.007297       0.003830      0.003851      0.004830      0.004919  \n",
      "1          0.007797       0.004922      0.003119      0.004019      0.005827  \n",
      "2               NaN            NaN           NaN           NaN           NaN  \n",
      "3               NaN            NaN           NaN           NaN           NaN  \n",
      "4               NaN            NaN           NaN           NaN           NaN  \n",
      "...             ...            ...           ...           ...           ...  \n",
      "1188            NaN            NaN           NaN           NaN           NaN  \n",
      "1189            NaN            NaN           NaN           NaN           NaN  \n",
      "1190            NaN            NaN           NaN           NaN           NaN  \n",
      "1191            NaN            NaN           NaN           NaN           NaN  \n",
      "1192            NaN            NaN           NaN           NaN           NaN  \n",
      "\n",
      "[1193 rows x 31 columns]\n"
     ]
    }
   ],
   "source": [
    "rois = []\n",
    "df = pd.read_csv(tbpath + 'fit_table_reweighted.csv') \n",
    "rois += list(df.roi.unique())\n",
    "  \n",
    "rois = list(set(rois))\n",
    "rois.remove('US')\n",
    "rois.remove('Gambia')\n",
    "#sort within US and among other coutries then union back\n",
    "roi_us = np.sort([i for i in rois if i[:2]=='US'])[::-1]\n",
    "roi_other = np.sort([i for i in rois if i[:2]!='US'])[::-1]\n",
    "rois = list(roi_us) + list(roi_other)\n",
    "\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AA_Global\n",
      "[0.01823967 0.04241168 0.06967841 0.11933414 0.37706928]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAIMCAYAAABc9G6eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7gkVXnv8d8LagAJoKImooJCFJN4QbzFEzODYpRElHg55II6KhIFfc7g7QFOOLNHkpiLTsgTRdneUMEcjok3okRFZtAcJGZUojEoeghbvAaMg8hFRdb5o6qhp6e616qqVWtVdX0/z7Oehu6q6tpr9u5e9a633mXOOQEAAHRlt9wnAAAAlhuDDQAA0CkGGwAAoFMMNgAAQKcYbAAAgE4x2AAAAJ1isAEAADrFYAMAAHSKwQaQkRWebWbnm9nVZnZT2a4qn3uWme2+YP+9zexHZubMbLvnva4ut5tuN5nZ18zsbWZ2aPyf8Pb3PrJ8v1d39R6pTPou4vHuZGavMrMvlv8e3zezj5rZb8zZ/h5m9kMz+2CscwC6ZlQQBfIws/tKer+kR0tykr4o6euSbpP0QEmPlGSStjvnHj3nGC+Q9I6ppx7qnPu3OdteLelASR+T9N3y6XtLeoyku0u6RdJTnXOXtPrBdn3f3SVdLukekg52zt0c8/ipTQYazjmLcKw7S/qopCMlfV/SVhX9tK7c5AXOuXdX7HeGpD+S9ETn3Na25wF0jcEGkIGZ7S/pc5LuL+liSS91zl05s819JJ0m6fecc/eYc5xLJP2GpG9Luo+kNzjnXjVn26tVDDaOcM5tm3r+bpI+JOkJkr7qnIsa4ZgaEL3SObcl5rFziDzYOE3Sn0j6gqQnOed+UD5/pKQLJf1M0oOcc9+Y2W8/Sd+SdIVz7lFtzwPoGtMoQB5vVjHQ+JSKaMKVsxs4577tnHuZpGOqDmBmD1QxQLhR0gvKp//AzO5U50TKL7jXlP/74PK4Mb1c0k8lvSfycQet/Hd6Rfm/J04GGpLknLtI0tsl/Zyk/zG7r3Nuh6QPSjrczH4twekCrTDYABIzs1+S9Kzyf09yzv100fbOuU/Peen5KqZZ3u+c+7ikKyX9gqSnNjit6amXezfYv5KZPUbSYZI+5py7dur5O5V5BzeZ2V1m9nn5VE7J/WZeO6p8/m8r3mtvMzvNzD5vZjeUx768zIe4y+z2TfdZ8LPe1cwuLM/vA2a2p2eXx6uYMrnaOXdZxev/u3x8xpz9J9MrL61znkAODDaA9J6mYpDwr/PyK3zMzFQMNiTpnPLxXeXjhgaH3Gfqv/+zyTnNMfmivHj6SefcrZI+LWlPSbNX5k+a+u8nznltp+OVg5LtKqYkflFFxOhiSfeV9JeSPlYxqKm9zzxmdi9J21QM9N4s6VkBuSmHlY+fm/P65PmDzeznK17/tKRbJT3NzPgsR6/xCwqk98jyceHdIx5HqMi/+IaKpEKpuNK9TdLRZnb3msf7rfLxK5KuanFes9aXj1VX7pMBw+0DivJLc52kL6tImn3SzD6TbS+e2sckvU/SgyW9QdJBzrnfds49TdIhKhJi16tIqGy8zzxmdoikSyU9StLpzrkTnXO3+fZT8e8nSddUveicu0HSD2e2nX79JklfknQ3SQ8LeD8gGwYbQHr7l4/XLtxqsQ3l43tcmeXtnPumpE9Kuouk3w85iJndy8yer+JK/gZJx0+OV7HtIWa2YmZ71TjPR5SPV1S8tstgQ8VAbD8VA4EvaueByN0lPVzSN5xz/29qn6MkPVbSJZJe7Zz78eSFMrfhBZJ+IunEcpDRdJ9dmNmjVQw0DpT0QufcH8/btsLe5eONC7b5UflYFdmQ7ujXR8x5HegFBhvAwJQh9UnOx7tmXg6ZStk6VSvieyqmYW6S9Ejn3P9dsN9TJG0sr6hDzvOukvZScUfF9RWbXC7pvyQ9ttxW2jlycbGkA8zsweVz61V8Zu00haJi4CBJf1c1UHLOfUfS11TkR/xSi31mf77fUhFV2kvS051z76zarmPfLx/vleG9gWAMNoD0risf79lw//+u4gvuUufc12Zee7+K0PvhZvarc/b/mIpByXtU5Bn8TMVts+81sz0WvO/DVITtQ+1XPv5ozhe6K9//ziruqpGKaZObVEy7XDz13PTj7GBjcvfM31QULZsMqn6l3OaeLfaZ9SFJd5X0bOfchXO2WWQStbjrgm0m0Y8b5rw+mWbZb87rQC/UukUOQBSfl3Scijn+JjaUjwea2T9VvO6mtququfFnM3U2HqLiCv3Rkv5Ud9yOqaltvqPiTpfb60xIeoZz7sMLznNH+bi3mdmc6ZmLJT1T0hPN7JOSfl3Sp51zPy1riNyqYpBxliryNUq7Tz1fmf8wZRIJaLLPrHNV9PFfmNnnnXN1E2vXysf7Vb1YRrD2mdl21r7l4445rwO9wGADSO8jKpISH25mv+Kc+3LojmZ2sIovZEk6oGzz/IGZnVLe+TGXc+4KM3uxpA9LepmZvdE5N5sk+jsq8kHepCJ6InmiHM65G83sRhVX7vuq+gtxOm/jcSoiNheX+99gRQn29WZ2gKRDJV3pnPvWzDEmg4X3OufevuicWu4z64Uq6oe8WNI2M3uic+67nn2mfaF8PHzO65PnryqTRatMEoFj3kEERMc0CpBYWcDrA+X/vqksWT2Xmf361P9uKB8/4JyzeU1FvkFwzQ3n3AUqoht3VvUdGN9UMRD4iHPusrItSmycuLx8/OU573uFpO+ouA10kocyHbm4WMUX6isrXpv4x/Lx2QHn02afnZSRmj+U9EZJD5F0STkoCnWpiqjJQWb2uIrXf7d8XLQGyqRfv7BgGyA7BhtAHi9V8QW+TtKFZaGvnZjZvc3sTBW5AZPbQp9Xvnyu5/jnlY8bapzT6eXjc83sATOvPbx8/GKN40lFToZURC3m2aris+gPJf1AxTTTxCfLx0nhqqrBxgdUfNk+1cz+ysz2md3AzA4ys+Na7rMLV3i5pC2SHiTpU2Z2/0X7TO17a7mfVAw6b8+7KMuVv0jSjyX9ddX+5V1BD1WRZFv33wVIyzlHo9EyNBXlyreryLG4TcWX3/sknS/psyoSN52ky8rtjyz//weSfs5z7EPKbX8s6e7lc1eXz61fsN+F5TZvm3n+NElrDX7Gx5TH+4cF27yo3MapiNhMv7aHpJun+mj/BX05qc2xQ8UtreepGKhdOd2PbfYp93MqAxszz/9p+drVkh4Q2D93lvSJcr/ryn//T5b/9rdJev6CfY8q93t37t9lGs3XiGwAmbhica3HSDpW0t+ruM3yaZKerqIWx/tU5Eo8vtxlQ/n4PjdVF2LOsb+uYsASXHOjNJlCed5MdOMRkv61xnEm5/FZFYOop5RVNqvMTptM73+LiukGSfqSc+46VSj78lGSNqoovf4wSc8pn/u+pD+WdELbfRZxzp0mabOKmhufqopWVezzUxUF1V6jYiXep6moNfJxFYPC2Vubpz23fHxz6DkCubDqKwAvM7tc0kVuzoqynn0nq76+2jn3+ugnN0LGqq8YGCIbAELskPQEM1tnZo8zs929e9zh3SoiB68MWJwMYV6pImH31blPBAjBYANAiNNU3MJ6kYoVXH8WumO57StU3B3zsm5ObzzM7B4qlp3/kHNuq297oA+YRgEAAJ0isgEAADrFYAMAAHSKcuV+zDMBAMbEYh+QyAYAAOgUg42eOvXUU3XqqafmPg0AAFpjGqWnPvOZz+Q+BQAAoiCyAQAAOkWdjQpmdoLKNRHOPvvsw084oXp5BLOwHJqNGzdqv/32829Y2uS2aLO9Inj7unbs2KEzzzyzs+MDQ8fnIkYueoIo0ygVnHOrklYn/+vZ1nu8M888U9dff32EM4vjwAMP5MMUmCP0IgJAOCIbfq07yMxqf7kfd9xxkqRzzz134XZdD2S6jrLEsO+++2rjxo25TwNLrMnfMDBgRDZSmJlG0bxplJrHbLTfeeedt/D1ulM0y2htbY2rUQDoMSIbfnQQAGBMiGykUCey0eSKuk00IvW0BsmkGDMuxoA4GGxUqJMgWm5f6/h9SxhdhGRSjBVTc0A8TKP4ReuguoOMLqMYJFUCAOZgGiWFrqZRVlZWtGnTpqBt169fr03rFbx9XZs3b+bKDZKYKgDQPSIbfkQ2AABjEv1KlMGGn7eDurgHfxJNWV1dHVSOR4gh1O7AsDB4BqJiGiWFJnU2upqSeOtb30otDcCDWitAPF0EIYhs+NFBAIAxIbKRQkhkw3cVFRKNiDGdQB2M+BiAA0BcDDYqhNbZWPSllCrPgjoYcRGKB4D4mEbxy9JB69evlyRt27Zt4XYpk0f7nthJkiAARME0SgoxF2Jre6UcY7pmLEgSxAQXUUC/ENnw63VkAwCAyKizkcJMZOPwrhJEJ1JMT5BIihB8HgAQ0yhpDClBNBSJpPBhCgpAVxhsNNT2i3s6mtEmqkFSJGJhMAqgKww2KsSos+FddG1ly8LXjzvuOEnSueeeu/B9WFAN0xgwAOgjcjb8GnWQbxolVp4GkQ0AQGQkiGaQbCG21Hkefa+bMSQM+gAsERJEU8i1EBs1M4aLGh8AlgULseVBnQ0AwJgwjZJCSJ2NqW29xwuNWMSc1qCuRv/xtwegp5hGSSG0zsbU9gtfz1Fzg7oa/caUC4AxIbLhV6uDUtyFQjIiAKBD0a+Gdot9wGVgZieY2XYz2766ulr1+ty2trbW+flNkhFpNDMiJAD6j8iGHwmiAIAxIWcjBWu4xHwXV5lmlvyW2GWqv0GiLDAeXDz3F4ONCnUTRGf2jXIO05GNvi3qNiQkygLjwJRivzGN4le7g0gSBQAMGNMoKYRMoywaRbddhE2SzjjjDEnS6aefXvk6C7CNFxcIAIaGyIYfkQ0AwJhQQTSD2h1k1m5hthg5GqmTPBn8AHdo+xkAZMY0SgpN70aZOUbj9x/igmwsRAYAmIfIhh91NgAAY8I0SgoWeSG2uvoe2VimOhzIg/onQH8555hGSSH2Qmx1UVcDy476J8C4ENnw66yDhpgIOovEUABYOkQ2UqibINp0KmVRPY7Xv/71kqRXvepVnoP4a3Z0iXofQHNc7GEsiGz4EdlYgMgGACwdEkRTCEkQ9V3NhyR5th0okGSXH38/AJYQ0ygphCaILvqiSZHkSZJdXkwfAUAYIht+yToo5gAl9/TKBNMsADA4RDZSaFpBtO2V7nTCaOuiXpkTRydIIMUy4KIMaIfIhh+RjRaIbADA4JAgmkFwB8VcfGk2spGy0FdfBiqhGNAgFxZcw5JiGiWFNguxxZ4ymByv7yXMc2IROADoNyIbfnQQAGBMiGyk0HUF0aZRiq6nN6jbMQxcIAAYGgYbFeouxFbuE3z8vi60Rt2O/mO6CMAQMY3i12kHtRl4pE7kJBETAEaBaZQUmiSIxl6MLajORuJaGtTMGB8uRgDEQGTDj8hGicgGAIwCdTZSiLEQm1Q/EZQEUGC8+CxGjzCNkkKMhdik/iWCkgAK9BPTk1h2RDb8Oumg1AORoVUFHRumqAD0SPzRr3OONtNUTKFsl7T97LPPdvOoGIg0aisrK3OP65xz1157rbv22msXblPLpn3iHQvRraystPp9otFotFjNdfC9SmTDr5MOIrKBaUQ2APQIORt91GQxpnlfLH3L80iJL1wMDQuxAWF2y30CfWRmJ5jZdjPbvrq66t+h2CdKW1tb6/in66/Jgmo02lAagDBMo/hl6aCgol4AAMRHnY0ULEKdDV+NjbY5FNTMiIe/AQDYCTkbKbgIdTa6zr2gZkYchMIBoHtENvw676A6A5OmERGSLwEAgYhspGANFmKb2b/W9lWLsV1zzTWSpPvd734zGzdbfG0oi6gx+AWA5UNkw4/IBgBgTEgQTSEkQdSzf9TzqbugW4g+Ffki2XXc+AwCeodplBRCE0Q9x2h1DtO3vi57oS+SXcdrCFN7ANojsuEXtYNiDRz6FJmoiykdAOg1FmJL0dThQmy+Bdgm1q1b59atWzd/gwEvrMaiYzQajRbeMmAhtgyidhCRDSIbANBzJIhmENRBZs0WZGo7+Mg16GDAACzW9DMB6AESRFNoWmejSbJbF3eapDBZNA0AAB8iG35ZOoiF2AAAmTCNkkLdOhs5Iho5pk+oh4G6+HwBBolplBRcgzobdT9Uh1g7g3oYqINpNgATRDb8anWQbxBBuXEAQM8R2UjBlyC66IqtalG1nTcIW0htNmdjKAupLRMG4gAQB5ENPyIbAIAxIUE0g7kdFHIffWhuRuyETwYn6BI1JIClFn2wsVvsAy4DMzvBzLab2fbV1VXftgvb2tpaorPe2aQOBo3WRQOAOohs+PWqzkaqu1iGUg6dCA4AREeCaArWsILo1P4xz2Wn//cmoMYSmMiaG4mzyIGLNKAeIht+RDZ6jMgGAERHgmgGjRNEUyeH8sU7bL7fJwBIhATRFCxSgmjq5FCSQofdAGBZEdnwy9JBl156qSTp8Y9/fI63BwCMF9MoKVjAQmyhV6J1F1zrIleCBdT6j79DAD3C3SgphC7EFvIF0YcF11hArd+YQgGw7Ihs+HXSQbEHIUO5e2SCZFYA6C0iGylYjTobTa9KffUy5t36Ov+Aw6iLMUF9DKB7XEyiL4hs+BHZ6ACRDQDoLRJEM2i1ENu0NgOMXIMJBgUYmrp/lwB2QZ2NFCziQmzTLdeibG1Qu4M2tAagf4hs+NXqoJDoRYwoBREHAEBHSBBNwTwJoouunoIWSouQzEmCZT4M0AGgHiIbfkQ2AABjQoJoBjt1kNn85LNFA422AwwGF8O06PcFAHqKBNEUzJMgOi8xrcsEUBI1h9kAAEQ2QmTpoNpFvQAAiINplBQsYCG2gGO0Ooe6C7gt0oeCXywGh7HiMxYDxN0oKYQuxBZwnMbncMoppzTet49YDA5jxFQaUCCy4Zesg2KVMO9DJIOEVgAYLCIbKfjqbMzZp/X7BtXoCDpQ/kXZqAOCPuMiC0iLyIYfkY0GiGwAwGCRIJrBwg6qU0ehyWCiy4EDA4Jho4YHgI5QZyMFX52Niu2DWt8WYqN2x7AbAAwFkQ0/6mwAAMaEaZQULLDORujVZd2aGalyLqh9AaAtvkOWEnejpBBaZ6PLXI0UqH0BAEiByIbf3A7qOjlUShPlIFEUbZCoCiwdEkRTqJMgOtTk0GkkitLaNADwIbLhl6WDzj//fEnSsccem+PtAQDjRYJoChaQINr2ii7mQmtALH0oCAcgr02bNpEgmkLsBNEqfU0aBQAgNiIbfq07qM3Aoi9XmiSRAsBoENlIwQIWYqszjdJkgbVJUa9N65V9UTWJhdUAjAsX4nER2fAjsiEiGwAwIiSIZhDUQWbNag20zd3IMRhh4IGxa/r3DgwE0ygphEyjzNmv9nsN8a6USV0OAABCENnwy9JB55xzjiRpw4YNOd4eADBeTKOkQJ2N4ehLTgsALAvqbCQSUmej7SCNOhsAgLEgsuHXaYLotNgDkD5e9ZNcijEioRQDQ2QjhZQJotPGMLVCcikAjA+RDb8sHTQp6rVt27Ycbw8AGC8SRFMISRAttws6Xp2IRVdTHzt27NCZZ54Z/bjoF/6eAUTANEoKoQuxldt6j9eHZNADDzyQL6Ilx/QUgL4isuEXpYOaDjhSJXmSuAkAKBHZSCE0QbTLxdhSL8TGQmtAOlzkYWyIbPgR2QAAjAkJohks7KC6988PdQVYBiMYEupaAK1EH2zsFvuAy8DMTjCz7Wa2fXV1NWT74La2tpbgJ4hvUh+DRhtCA9AvRDb8qLMBABgTplFSsMA6G1PbR33/IVQS7WMpdORHPRdg+Jxz3I2SQp06G1P7RHv/PtTlAJqgnguAKgw2EGQ2kkFUY5hI9AWQA4ONCtZgIbaYUylDmEbBMLEQHgCfLqKT5Gz4kSAKABgTEkRTCE0QrXOF2CRakTIJk8Q+pMDnDTAIJIimEHshNqn/SZ8k9qFrTN8A40Vkwy9qBzUZdKSIcJA4CAAoEdlIoU6CaN2rtboLshU7bel8MTYWYsM8XJAAaIvIhh+RDQDAmJAgmkFwB5m1X/wpVm5H7gqfDF4wT4y/EwCdYiG2FKzmQmwz+7ZqQ12obRYLt9HmNQDjQ2TDjzobAIAxIUE0BWtQQdRzvKj7dl1hNPcUTFPUCsFQcJGHsWGwUaHJQmwBx6y1/aLIRt9rduRCrRAMAVNJGCOmUfx26SCzxQluvsFArMgBSZjp+P7NAWCJkCCaggUkiC5KgEuV5EkSJkmNADAERDb86CAAwJiQIJpCnQTR0Kve0KTOWFMsJEs2xwAcAOJisFGhboJoyJdT6qROkiWbYcoEAOJjGsWv0w5qOwjJfZsqSaoAsHTiX3U552gzTcUUynZJ288++2wXQsWgpHZbWVmpPN66devcunXr/G+8aZ+g8+vKyspK45+dRqPRfA1ZRP9eJbLh12kHEdkAAPQMC7Fl4O2gJjUYhjbIYFCB3Kh1AiRDnY0UQupsVOxTqw1twTVqetByNwDDRWTDjw4CAIwJdTZSsMA6G76rrZDaGtTVGAcG9QDGjMFGBVejzsaiL5GUtTWoq9FfTAEAGDumUfxu7yCz6gS1OoOK3HePjBlJrgAQhGmUFBZNo1RdpYaWIkdekyRXAMB8XQQhiGz4Zemg9evXS5K2bduW4+0BAONFnY0UZiIbh1cliMZIDpXiT6uQKIou8XkBjALTKCmEJoj2JTl0Gomi6ApTUACaYrDRkO8LfePGjckGHCQ+IgUGsQCaYrBRIVadjZWVFW3atGnxm61sqdxm8p4hFUw3b97MVecI8eUPYCjI2fBr3EEhkY0YORtENgAAEZEgmkHtDppXj6OOyd0oxxxzTPSpmL7U+mCQhBh/KwCiI0E0hdBpFM8xopzLYYcdtrQ1PKh7AQDjQGTDjzobAIAxYRolhZA6Gwv2bfXesaqR5pwqodYH+obPOaAWplFSqLMQ25z9G793rvocMVHrA33CVB2QH5ENv2gd1GYgkTpSQfImAIxW9BH6brEPuAzM7AQz225m2xfVuTCzWm1tbS3hT9HOJHmTRsvZACwHIht+JIgCAMaEBNEUbEGC6LyrrUWJnW2nQEi4BIB2+K6rhQTRFHwJolW/tF0mdpJwCQDNMSWXH5ENv1YdVGcQ0mUSKAmfAIBARDZSsEgLsUmBi7HdvvEdi7KdeuqpkqTXve51Yft6sFgb+oCLG2CciGz4EdkAAIwJCaIZeDvILGwxqSZ5HbkXTWOQAgCjwzRKCqHTKDP7eLeJVYo8JRZLA4Bx6SIIQWTDjzobAIAxYRolhUV1Njz7tX7vmNGP3FMws6gXgjHhsxUDxjRKCm0WYmv7AbMMC7HNQ70QjAVTj8DOiGz4ddJBMQYVuSMXJI8CwFIispFCnQTRplcwvvobZ5xxhiTp9NNPX3CQLeE1PDpA7Q7kxIUSMBxENvyIbMxBZAMAlhIJoimEJojWvapvmvyZalBBAie6wucMMChMo6RQJ0G0zodo35M/SeBEF5hqA0Bkwy95B6UalOSehomF6RwAiIrIRgp1K4jGvnJbWVnR1q1bJXVc1CtzgmksJKoCy4ML4OVEZMOPyEbPEdkAgKhIEM1gbgeFLsBWpesBxbIMJOpg0AEAUTCNkkKKOhtDXJSt71g0DgDaYyG2POggAMCYENlIoasl5ruw7BGSMU4HjQm1XYD+6SIIwWCjQtOF2GL+A4UuMd/32h3AItR2AcaBaRS/aB3UZmCQ4wqfhEsAGCWmUVLoKkHUt/jatF0iGxlqYlC/AiG4YAHgQ2TDj8gGAGBMqLORQuhCbFPb136PNomdqQceJPGNE58NwGgxjZJCkwTRuh/MQ0rsJIlvfJg+AxATkQ2/oA6qGjwsikAwRQEA6CkiGyn4EkSrrvoqkz8XJHUuc/IlA1gAwDQiG35ZIhuhdTYAAIiMBNEMWndQkwXb6gw2WCW2GlNVWEZtFoAEAjGNkkKTcuUBx+xsv2UvWd4UC7MBQD8Q2fDL0kFMowAAMmEaJYUUdTYmYkUlck9xUIsDY8VnKJYQ0ygppKizMTGkehuLUIsDY8Q0HRCGyIZftg5KMRDJHREJQaInACRFZCOFpgmisa9y6izc1vxN0i/wVtcy1yTBcHBhBjRHZMMvW4Lo4x73OO25556dvg+RDQDADBJEUwhJEF10pe1L+ozxBU9C5vLjbxNAJkyjpBCaIDrvyyBFrgUJmcuNaSMAy4TIhl/yDoo1WMk9RcL0BwAMEpGNFNpUEI1xRbqysqKtW7dKalnUK3PyJ4mdSI2LJ6CfiGz4EdloiMgGAAwSCaIZBHVQ28WRYud55B5oLMIgBEAdLD6XHNMoKeSoszGmxdRYIA0AxoXIhh8LsQEAxoTIRgptl5iPedU+e6xUEZA+T8PMQ+0RDA0XexgLBhsVmizEVnGMVucwL7KxLAu3dYHaIxgSphIxJkyj+HXSQU0GDV1GG0jaBACUmEZJock0St2rFN8ia5WRjQ7rZlATo5+4GACwDIhs+BHZAACMCXU2MpjbQb57v0MGFG0HEAwS4uOefgAjF32wsVvsAy4DMzvBzLab2fbV1VXftnPb2tpa5+c6qVlBi9cAAHER2fBb2EGLohdtohY333yzLrvsMupsAABSYxolhZkE0cMnCaJVV70LEz1X9pVWmt2munnzZq2srDTaF3HxNwJgZBhsZEBkAwAwJgw2MmicIDpPjMJcfajwSXIqciOZF+gECaIpxEoQzZk4mgLJqbTcDcAwENnwy9JBF110kSTpyCOPzPH2AIDxYholBZuTIBq4b+v3j7XYWu6pFhZGQx18FgG9QbnyFNouxNb2Q3NZFltjYTSEYkoEWG5ENvyidVDfSpRLJHkCAHZBZCMFq7EQW50rMt/ia9MmC7FtWq/OFl+TWIAN7XCxAiAEkQ0/IhsAgDHh1te+6TIqwEADIYhMAeg7BhsVrEadjXL7oNa3+hrUyViOBgB9xzSKX6MOSrG8/AQRCgBARCSIpmABCaK+K8qgZNCVLVGSP0nyRAguLADkQmTDj8gGAGBMqCCawS4dZDZ/8SffIKPpAIMBRbcW/ZsCwMhwN0oKFg6mC9MAABo+SURBVJAgmjoJlGROkiwBYKiIbPhl6aBJUa9t27bleHsAwHgxjZKCNViIre3VcYzF13ItvMaCa4iNzyUgK+5GSaHpQmxtPiCHvPgaC64hJqa1gOVDZMMvSgf1sVT5BMmnAIApRDZSCKmzUW4XfMw6i7DdsVOcOhw+1OnA2HCRBaRFZMOPyAYAYEyIbPTNEOszMLhAE0P8XQfQD9TZqBBSZ2Nme2/r0yJs1OygNWkA0BTTKH6DrrOR6y6XXLfhtkHEBwAkMY2ShgUmiLY4fifbVmmUmBpDouTWmEiUxQQXYUBcRDb8iGw0QGQDAAaLCqIZ1Oogs3ZJdDEHB7m/8PnyRlNt/44AtMJCbClYzQTRiv0btz4lkrZFIiqtaQOwXIhs+GXpoAsuuECSdPTRR+d4ewDAeDGNkoIFLsTmuwILWVwtxlQHC6H1F39fAAaIu1FSqLMQ26Ivk1TJmSyE1k9MBwBAgciGX6/uRkkxgMmdWLoISacA0DkiGylYhDobsa5qZ4+TpG5Gj2tkUAsDXeCiC+gWkQ0/Ihs9QmQDADpHgmgGlR1kNr8OgG9A0ObLnC/bZhb9ewEAdkKdjRQssM7GvBoBXdbKoHYFdRsAYGiIbPhl6aDzzz9fknTsscfmeHsAwHgxjZLCojob866SfTU12uZBUEsDwCJ8liMi7kZJYVGdjaZ5Gm1RSwMAMFRENvyCOqhtAmKswUof7iQhiRXLjoRjLDkSRFMITRCt2K9xYwE2Gm04DUA9RDb8elVnAwCAjpEgmoJFWohNCluMTYo7/UEyaX/w9wVggEgQTSHWQmxSusXYppFM2g+E2wGgQGTDr1EHhQwy2kYzSMQEAHSAyEYKFrAQm++qNWjBtAULnoXkbLAoWbcYiANAHEQ2/IhsAADGhATRDGp3kFn7e/AnkY1jjjkmes5HH2pxzMMACssoxmcCkBDTKCmETKMEHCPKuRx22GFBd7Msi0mNDgDA8iCy4UedDQDAmDCNkkKOOhtSt9Mb1N5AE3w+AKPENEoKQ6+zUYXaG6iL6SwAsRDZ8GvVQU0Tw8466yxJ0oknnrhwu5SDmT4nli5C0imGgCRS9AiRjRRiJIjOHK/xvieddNLC1+tM04wVSacAkBeRDT8SRAEAY0KCaAqhCaJz9o1yDl1GLPo4HUICKyQSUoGeYBolhToJonP2b30Op5xySutjDAkJrGCqC1heRDb8GndQk+TNlFEHEicBABWIbKQQYyE2KXAxtl12mr84W2ws5IYh4gIJGB4iG35ENgAAY0KCaAYLO6juvfFt62LkSO5kUIK6qBkBDFr0wcZusQ+4DMzsBDPbbmbbV1dXQ7YPbmtrawl+grgmdSpotNAGANOIbPhRZwMAMCZMo6RgkRZiC6mV0XZahPoUfvyOA0At3I2SQqyF2FKsW0J9isUI6QNAfkQ2/DrrILP5SXRHHXWUJOnCCy/s6u1vl2tl2j5WMg1BwizaWvS3D/QAkY0ULPJCbJ73avV6DCzmVg8LuwFAPUQ2/OggAMCYENlIoU5ko+0VbqqowlCnLIaAJF0Ay6SLIASDjQp1F2Jr8w+TK18C8ZCkCwCLMY3iF72Dmg4wUkQnSH4EgNFjGiWFrqdRQhZoqyzqlWCRNhZnGzcuPgB0gciGH5ENAMCYUEE0g7kd5LtXvs6gIvZAgkHDuFC3AUBELMSWgtVYiG3RYlQ5F11j8bRxNQDoMyIbfnQQAGBMSBBNwRLW2ZjV12qe1OlAHdQeAYaLOhuJpKyzMYu6G1gG1B4BMI1pFL/KDjKbn5AXMmBoEykg+bO/Fv1eAMBAkCCaggUmiM5L1us6MZTkz/42AMCuiGz4ZemgyqJeAAB0jzobKcwkiB5elSAachVbJ9mz6wRMEvaAMHwmAtyNkkRogqjvQ6lPyZ4k7AF+TIUB3SCy4Verg+YNMEIjFyR/AgAyI7KRQkidjXlXQHMXWQtcRG2MC6Ex4AWA5UZkwy9LZIMEUQBAJiSIZhDUQW3qK8TK7ehTlU+mg4DlQf2Y0WEaJYU65cpn9mv0fn0tUd7GpBYIAABENvyoswEAGBMiGymkjmyEHjN1BKRP0zKpUI8EbXEBB+yKwUaFuguxzewb/4RKfarbsayoR4I2mDoEqjGN4he1g9oMGHJFGkj2BIBRYSG2FEIXYiu3rdW6XqStCyz8RkvdACwXIht+JIgCAMaEOhspWISF2EKSOdtOi5DMWB+/7wDgxd0oKcRYiC1FMifJjPUQngeAPIhs+EXroLoDkK4TQkn8BABUILKRgtWos1HnannuIm0VNm7cKO2n4O2bGOOib8uICwYAfUdkw4/IBgBgTEgQzWCXDjKbvyiRb0ARYwDBIAFIZ9HfO7CkmEZJIWQaZd70Q4qS4ixyBgAYEiIbftTZAACMCZGNFBZFNupGNJpMm+zYsUOXXHLJwvfDuHBRAGDIGGxU8NXZqPrgj1lX48ADD9S6deskEdkAA04Aw8c0il/WaZRjjjmms+JgfV5CniRYAMiGaZQU6tTZqNg36rkcccQR3dXaWNnSaR2PNqgBgti4sALyIbLhl7WDuix7TmQDAFCBOhsphCzEVrFP7fdpc5tsroECi78BcfEZjB5iGiWF0IXYKvar9T4pFmuLjcXfgHiYKsRYENnwG2Wdjb4OhPo89dMWU0cAeoLIRgptEkRrvEfU7WKrs2hcUj1Oam2LpFh0gQtK9AGRDT8iGz1CZAMAOkeCaAbeDmqyUFOML/O+fPHyJQkMC4vLwYNplBSaTKPUDX+nWLAtFRaGAwAsQmTDjw4CAIwJkY0UUkQ26uhrFKQv0zhjRc0TAF3oIgjBYKNCqjobixx11FGSpAsvvLC3yZrIi5onAIaCaRS/LB10yimnaM8998zx1rUR4cAiJBADg8M0Sgp1p1G6mEJZWVnR1q1bh7HE/BLXvkB71A8BhqWLIASRDT8iGx5ENrAIkQ1gcKizkULdhdhS3faa4kudpMPlxN85gBqYRkmhSYJonQ/zPid8knS4fJjCAJAbkQ2/Vh0UOrCIEbUgXA0AiIDIRgqhCaIhV4zBC5pFSLIkES8uBuIAEAeRDT8iGwCAMSFBNIO5HbRoMaOQQUbbAQaDi+6xYBWAEYo+2Ngt9gGXgZmdYGbbzWz76uqqb9vKtra21vl5ThZAo3XXAADtEdnwy9JB69evl6RhFPUCACwTplFSCK2zsejKN6SWRttpFGpixMPfAQDcjrtRUqhTZ6NNzkZb1MSIg+kSAOgWkQ2/xh3EnSgAgAEispFCrDobbWps3HTTTZKkvfbaK+CMqbHRBwzcAaAakQ0/IhsAgDEhQTSFOguxNYkmNF2IbVrOlVZJTEUMfPYAvcU0Sgp1F2Kr+6HZ54XYQpCYiraY8gPGhciGX+sOSjmdMo2pFQBAA0Q2Uoi5EJvULFE0RlEvkkaHgQE/gGVHZMOPyAYAYExIEM2glwuxMZDoBxZqA7CEWIgthSEsxMYibP1oAAA/Iht+LMQGABgTplFSCKmz4buq9dXSiJGfQb2L/uHvCcAS4G6UFELrbCz6YmEhtvFhWgUAqpGz0VDuL3kSRPsn9+8EAPQVkY0KMepseGtrVCy+Ns2Xs0ENjX5ggAEAfuRs+DXqIN80StucDSIbAICOkCCaQVAH1am3UCefo4sF1xio9AM1OgD0FHU2UqhTZ2Nmv6DWdQ0OH2p09KMBwFgQ2fBr1UEpS5UTsQAARMCtrymEJoiW2y48VpNF2KRmRb1IGu0GA3IAaIfIhh+RDQDAmJAgmsHCDgpN8mPQMSwkbwIYMRJEU6ibINq3pFASQEneBIA+IbLhl6WDLrjgAknS0UcfnePtAQDjxTRKChawEFu53cLj+BZjm4hZS4PF2ZYHf5sAMuFulBRCF2Irt537WorF2GaxONtyYCoHwDIhsuGXpIPaDEy6qDJahcRTABgFIhsp1KmzUbFvo/ecrcdRq86GZ1G3WKjjsTy4yACQEpENPyIbJSIbADAKJIhmMLeDFtViCBk8sPLrMFBzA8DIUGcjhTp1NnLW1aCeBjU3AGAIyNmoEHo3iu9qt+u7UbjzBAAwBAw2KoQmiPquer2LsLVM7CRhszkGaQCQDjkbfo07yBfZIGcDANBD5Gz0SduowmZ7hTa5LY33v/7667V582YqhvYYkScAYLBRKUaCaKokUYlE0T43AADTKCGydFCtol4AAMRDnY0ULOFCbG3zNlh4rR1+/wFgF5QrT2FIC7Fx+2tzTHMAQBrkbLTQhy95kkSb68O/HwCMAZGNCsnqbEjRFlGj5kZeDFwAYD5yNvw6q7MhxVtEjZobAIBISBDNYGEHmYUv0tUkh6PLFV0ZoKCP6vxNAegERb1SqFNno9w+qKWquxGK+hy0PjYAy4fIhl/tDmJ5eQDAgHHrawoWkCC66AosRmKor6gXCaHLjwsBAMuCyIYfkQ0AwJiQIJrB3A4yW5zIxqCjv3z/dgAwYiSIpmADWoiNJE+SEAGg74hs+GXpoIsuukiSdOSRR+Z4ewDAeDGNkoIFLMTmuzr2LcLGAmzDx98OgCXF3SgphC7ElnMRNhZgy4upGAAIR2TDL1oHUUEUADAARDZSCKmzMbN90HGD6m+UJnU2Nq1XlIXaqlCrYzy4qACQE5ENPyIbAIAxIUE0g9YLsdUZZMQeXDCgwDzUGgEwB3U2UqhTZ6PcPnu9jXmow0Gb1wAgFSIbflk66IILLpAkHX300TneHgAwXkyjpGABdTbK7Rod31eDo60u8zyWETVLAOAOzjnuRkkhRp2NRbquwYF6qFkCAN0isuHXOkF0oukgI0WkgkTSYavzewgAHiSIpmCRE0QnLWeiqA+JpMNuANBnRDb8snTQpKjXtm3bcrw9AGC8SBBNITRBtNzWe7w6CaFdTJmQALlc+JsF0DESRFMITRCd2n7h67kTQkmAXB5MmQAYInI2WuJLHCnx+wZgiIhsVLDAhdhCrzLrLMCmlS3atGlT1JwNFlwbBwYiAPqKnA2/1h2Uc20UidtaAQC1kCCaAXejAADGhDobAABgWDobbJjZBjNzZrY+x/6eY59jZoR0AABIIDhB1Mz2kPRCSc+W9FBJ+0m6UdLXJF0s6Z3Oua90cZLLok7uxhFHHFHkb2zeHP08yOEAAKQUNNgwswdK+gdJD5F0iaS/kvQdSXtLeoSKQcirzOz+zrlvdXSug3f99deH35Ui3X5nSmxdDGAAAJjHO9gwsz0lfUTSwZKe6Zz7QMU2e0g6WZmSKZfROeecow25TwIAgAhCcjaOl3SopL+sGmhIknPuFufc65xz3/YdzMz2N7M3mdk1ZvaT8vFNZnaPObvcycxWzGzNzH5sZl80s9+tOO5vmtn5ZnaVmd1sZjvM7ONmti7gZ2ysq/oV55xzTifHxXKhfgqAIQiZRnl2+fi2tm9mZvtKulTSIZLeIenzkg6T9FJJTzSzxzjnbpjZ7c8l3VXSWeX/v0DS35rZHs65c6a22yDp7pLeLembkg5QMVD6pJkd4Zz7dNvzBwAA9YUMNn5V0g+dc/8x/aSZ7S7pbjPb3uicu3nBsV4j6ZckneScmwweZGaXS3pj+frpM/vsL+lhzrnry23fIumLkraY2flT7/di59yNM+f4FklflnSqpODBRmgF0antg467srISegqd44oYAJBKyDTKPpJ+WPH8QyRdO9NO8hzrd8rtVmeeP7t8/ncq9nnzZKAhSeV/v0XFQGf91PO3DzTMbO9yWuZnkv5Z0mM957UT59yqc+5RzrlH+QYa5fZBrU9Cz5nW7wYAQxAS2fihigHHrP+Q9OTyvx8u6fUBx3qApO3OuVunn3TO3WpmV0p6ZMU+V1Q89+/l4wMnT5jZwZL+RNJTVNyWu9NbBJwbAADoQEhk498k7WNmD5h+0jl3o3PuIufcRZI+18nZBTKzvSV9StJTJf21ijyTp6gYDF2sDkqvdu2jH/1o7lMAACCKkMHG35WPx0d4v6skPdjMdoqolP//oPL1WQ+peO6Xp44nSU+SdB9JJzvnVpxzf++c+3g5ELprhPOeq6tQ9l577dXJcbFcmEoBMAQhg423SfqKpFebWVVOhRQeOfigpHtq14HLi8vnq26tfWl5F0vxRsV/v0TSDhUFxqQiN2OX8zCz31TNfI2+OOuss/wbAQAwAN6cDefczWb22yoqiL7fzLZJ+rik76rI5ThU0rEqvvCv8RzuLyQ9R9KbzOyRkr6g4tbXF0n6avn6rOsk/bOZvbP8/xdIur+k451zN5XP/VN5Pm8ws4NU3Pr6CEnPlfQlFeXVs9t3331rV+/sqlw5AACpBC8xX1YSnV4bZV8Va6N8XUVexNudc1+d2n6DpHdKOsI5t23q+XtK2izp6ZLuLel7kj4kaZNz7rqK/Z8s6QkqBhn3lnSlpNc55947c34PUzFYeayKQdTnVNxG+yJJz3fO2dS258w+twBLzAMAxiR6nmPwYGPEGGwAAMYk+mCjsyXmAQAAJAYbAACgY0yjeJjZCc652YqniIg+ToN+7h593D36uHtd9DGRDT9/vXK0RR+nQT93jz7uHn3cveh9zGADAAB0isEGAADoFIMNP+YGu0cfp0E/d48+7h593L3ofUyCKAAA6BSRDQAA0CkGGwAAoFOjG2yY2W5mdrKZfcXMbjGza8zsDWYWtBR92/3HoE0fmdmDzOy1ZnaZmV1rZjeY2eVm9j/p4zvE/D00s73M7Cozc2b2xi7Od6hi9LOZ3d3MXm9mXy+Pca2ZbTWzJ3R57kMR4TN5bzM7zcy+VH5eXGdml5rZBjOLXnZ7iMzsVDN739Tf+dUNj/M8M/uCmd1sZt8zs7eV6535OedG1ST9tYr1Tt6vYmn7LZJ+qmIxud263n8MrU0fSfozSTdIOk/SyyW9RNL55fH+VdKeuX++PrSYv4eSXl/2uZP0xtw/W59ahM+LAyX9h6Rry9/tF0o6WcUik7+b++frQ2v5ebGbpE+rWHX8HSrqQ2yU9M/lMf8898/Xh1b2xfclfULSf0m6usExTi6Ps63s59dK+pGkL0u6q3f/3J2QuMN/RdJtkv5+5vmXl534+13uP4YWoY8fJWnfiuf/uNz/Zbl/xtwt5u+hpEdKulXSKxhsxO/n8ovwGkm/mPvn6WOL8Hnxa+V2fzXz/F0kXSVpR+6fsQ9N0gOn/vvf6g42JO2vYpX3z0rafer5o8v+P813jLFNo/yeitXszpx5/q2SbpJ0XMf7j0GrPnLObXfOXV/x0vnl46+2PsPhi/J7aGa7l/v8o4qrSuysVT+b2W9I+nVJf+Gc+46Z3dnM9urkTIer7e/yPuXjt6efdM79RNJ1Kr4gR885d1XLQxwjaS9Jf+Oc+9nUcS9QMajzfuaMbbDxaBWj6M9OP+mcu0XS5eXrXe4/Bl310X3Lx+81P7WlEauPT5Z0qKSXRT275dG2n3+rfPyGmV0g6WZJN5rZlWbGhUmhbR9/VtIOSa8xs+eY2f3N7FAze52kwyWtxD/lUZr8O3ym4rXLJB1qZnsvOsDYBhv3kXSdc+7HFa99S9L+ZnaXDvcfg+h9VF6Bn64i3P/e9qc4eK372MweIGmzpNc6566Of4pLoW0/P7h8fKuku0t6voqcjZ9Ieo+ZvSDmyQ5Uqz52zv1A0tNV5CH8H0lrkq6QdJKkZznn3hr/lEfpPuXjtype+5aK6NR9Kl673Z1in1HP7SWp6pdakm6Z2uYnHe0/Bl300Zkq5mZPc859tcW5LYsYffwWFeHPLRHPa9m07eefLx9vkHREGdqXmX1QRd//qZm9yzl3W6TzHaIYv8s/UpGH8GFJl6oY2J0k6b1m9gzn3CcineuYTab/qv6tbpnZptLYIhs3Sfq5Oa/tMbVNV/uPQdQ+MrMzVIT5V51zr2t5bsuiVR+XIfwnS3qpc+6nkc9tmbT9Xb65fPzbyUBDuv1q/MOSfkF3RD/Gqu3v8kNVDDA+4Zx7tXPuA865t6vIlfmupLeWkVG0M/k3qPq3CvpcH9tg49sqwnJVHXaAinDeohF02/3HIFofmdmKpD9ScZvgS6Kd4fA17uNyny2SPirpu2Z2iJkdouIWTUnat3xuvy5OfGDa/i5/s3z8bsVr3ykf79bi/JZB2z4+WcWX3fumn3TO3STpIyp+rw+Kc6qjNknAPaDitQNU3JHy7YrXbje2wca/qPiZHzP9pJntIekRkrZ3vP8YROmjcqCxSdK7JB3vyvusIKldH+8p6Z6SflvS16batvL148r/Pz7qGQ9T29/lSdLjfStemzz3n21OcAm07ePJl19V9OJOM49o7l/Kx1+reO1xkr7qnPvRogOMbbAxKQ61ceb5F6uYbzpv8oSZHWxmhzbdf8Ta9rHM7H+pGGi8R9ILRz6nXaVNH98o6TkV7cTy9X8s///DnZz5sLT9Xf6ginyN46Yz9c3sF1XcSnilc+7rXZz4gLTt438vHzdMP1lG5p4h6QeSxt7HtUzd0XPnqac/pGJa8GXT01JmdrSkByrkuy93sZEMxU3+RndUqzte0htUVKvbpqlqdZKuLrqn2f5jbm36WEVil1ORVf48FVfa0+3JuX++PrS2v8cVxztIFPWK3s8qKi06FQmMr5B0Svm7/RNJv5n75+tDa/l5caCKypi3qbg4eYmk01RUbXWSTsz98/WhSXquiinpP1JRPuAHU///3Jltt5V9d9DM868sn99a/l5vVpGce4Wkvb3nkLsTMnT67mWnfVVFZu23VMxh7z2z3bwPj6D9x9za9LGkc8pf6HltW+6frw+t7e9xxfEOEoONTvpZ0jNV1CK4UUWk4+OS/lvun60vLcJn8sEqplu/qWKQ8kNJn5L0zNw/W1/a1ADC+5k6b7BRvrZBxbIRt6iYAnyHpHuFnIOVBwAAAOjE2HI2AABAYgw2AABApxhsAACATjHYAAAAnWKwAQAAOsVgAwAAdIrBBgAA6BSDDQAA0CkGGwAAoFMMNgAAQKf+P4hTuwhHhtKpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "font = {'family' : 'sans-serif',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 18}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "def simpleaxis(ax):\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    \n",
    "\n",
    "def plot_table_data(col,ax):\n",
    "    _025 = 0\n",
    "    _25 = 0\n",
    "    _5 = 0\n",
    "    _75 = 0\n",
    "    _975 = 0\n",
    "    w = 0\n",
    "    for j, roi in enumerate(rois):\n",
    "        try:\n",
    "            boxes = [\n",
    "                {\n",
    "                'x': 0,\n",
    "                'label' : '',\n",
    "                'whislo': df.loc[(df.roi==roi)&(df['quantile']=='0.025'), col].values[0],    # Bottom whisker position\n",
    "                'q1'    : df.loc[(df.roi==roi)&(df['quantile']=='0.25'), col].values[0],    # First quartile (25th percentile)\n",
    "                'med'   : df.loc[(df.roi==roi)&(df['quantile']=='0.5'), col].values[0],    # Median         (50th percentile)\n",
    "                'q3'    : df.loc[(df.roi==roi)&(df['quantile']=='0.75'), col].values[0],     # Third quartile (75th percentile)\n",
    "                'whishi': df.loc[(df.roi==roi)&(df['quantile']=='0.975'), col].values[0],    # Top whisker position\n",
    "                'fliers': []        # Outliers\n",
    "                }\n",
    "            ]\n",
    "            sigma2 = (df.loc[(df.roi==roi)&(df['quantile']=='std'), col].values[0])**2\n",
    "            w += 1/sigma2\n",
    "            _025 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.025'), col].values[0]\n",
    "            _25 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.25'), col].values[0]\n",
    "            _5 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.5'), col].values[0]\n",
    "            _75 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.75'), col].values[0]\n",
    "            _975 += 1/sigma2*df.loc[(df.roi==roi)&(df['quantile']=='0.975'), col].values[0]\n",
    "            ax.bxp(boxes, widths = 0.5, positions=[j], showfliers=False, vert=False, patch_artist=True,\n",
    "                   boxprops=dict(facecolor='none',edgecolor='gray'))\n",
    "        except:\n",
    "            print(roi)\n",
    "    return _025,_25,_5,_75,_975,w,j\n",
    "\n",
    "\n",
    "def plot_global(col,ax,q):\n",
    "    boxes = [\n",
    "            {\n",
    "            'x': 0,\n",
    "            'fontsize' : 10,\n",
    "            'label' : 'Global',\n",
    "            'whislo': q[0],    # Bottom whisker position\n",
    "            'q1'    : q[1],   # First quartile (25th percentile)\n",
    "            'med'   : q[2],    # Median         (50th percentile)\n",
    "            'q3'    : q[3],     # Third quartile (75th percentile)\n",
    "            'whishi': q[4],    # Top whisker position\n",
    "            'fliers': []        # Outliers\n",
    "            }\n",
    "        ]\n",
    "    ax.bxp(boxes, widths = 4, \n",
    "            positions=[-4], showfliers=False, vert=False, patch_artist=True,\n",
    "               boxprops=dict(facecolor='none',edgecolor='gray'))\n",
    "    return \n",
    "\n",
    "# roi = rois[0]\n",
    "# col = theta\n",
    "# print(df.loc[(df.roi==roi)&(df['quantile']=='0.025'), col])\n",
    "\n",
    "# df['quantile']\n",
    "\n",
    "theta = 'car (week 0)'\n",
    "f,ax = plt.subplots(1,1,figsize=(8,8))\n",
    "\n",
    "_025,_25,_5,_75,_975,w,j = plot_table_data(theta,ax)\n",
    "q = np.array([_025,_25,_5,_75,_975])/w\n",
    "print(q)\n",
    "plot_global(theta,ax,q)\n",
    "plt.subplots_adjust(hspace=0.2,wspace=1)\n",
    "ax.axvline(q[2],color='k',linestyle='dashed')\n",
    "# ax.yaxis.grid(True)\n",
    "simpleaxis(ax)\n",
    "ax.set_ylim((-8,len(rois)))\n",
    "ax.set_xlim((0,1))\n",
    "plt.subplots_adjust(top=0.9,bottom=0.1)\n",
    "# ax.set_title(theta)\n",
    "ax.set_title(r'CAR$_t$ (week 0)')\n",
    "plt.savefig(productpath + theta + '_forestplot.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
